home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / scripts / control / dlyap.m < prev    next >
Text File  |  1996-07-15  |  2KB  |  105 lines
  1. ## Copyright (C) 1996 John W. Eaton
  2. ##
  3. ## This file is part of Octave.
  4. ##
  5. ## Octave is free software; you can redistribute it and/or modify it
  6. ## under the terms of the GNU General Public License as published by
  7. ## the Free Software Foundation; either version 2, or (at your option)
  8. ## any later version.
  9. ##
  10. ## Octave is distributed in the hope that it will be useful, but
  11. ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12. ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13. ## General Public License for more details.
  14. ##
  15. ## You should have received a copy of the GNU General Public License
  16. ## along with Octave; see the file COPYING.  If not, write to the Free
  17. ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
  18. ## 02111-1307, USA.
  19.  
  20. ## Usage: x = dlyap (a, b)
  21. ##
  22. ## Solve a x a' - x + b = 0 (discrete Lyapunov equation) for square
  23. ## matrices a and b.  If b is not square, then the function attempts
  24. ## to solve either
  25. ##
  26. ##  a x a' - x + b b' = 0
  27. ##
  28. ## or
  29. ##
  30. ##  a' x a - x + b' b = 0
  31. ##
  32. ## whichever is appropriate.  Uses Schur decomposition as in Kitagawa
  33. ## (1977).
  34.  
  35. ## Author: A. S. Hodel <scotte@eng.auburn.edu>
  36. ## Created: August 1993
  37. ## Adapted-By: jwe
  38.  
  39. function x = dlyap (a, b)
  40.  
  41.   if ((n = is_square (a)) == 0)
  42.     warning ("dlyap: a must be square");
  43.   endif
  44.  
  45.   if ((m = is_square (b)) == 0)
  46.     [n1, m] = size (b);
  47.     if (n1 == n)
  48.       b = b*b';
  49.       m = n1;
  50.     else
  51.       b = b'*b;
  52.       a = a';
  53.     endif
  54.   endif
  55.  
  56.   if (n != m)
  57.     warning ("dlyap: a,b not conformably dimensioned");
  58.   endif
  59.  
  60.   ## Solve the equation column by column.
  61.  
  62.   [u, s] = schur (a);
  63.   b = u'*b*u;
  64.  
  65.   j = n;
  66.   while (j > 0)
  67.     j1 = j;
  68.  
  69.     ## Check for Schur block.
  70.  
  71.     if (j == 1)
  72.       blksiz = 1;
  73.     elseif (s (j, j-1) != 0)
  74.       blksiz = 2;
  75.       j = j - 1;
  76.     else
  77.       blksiz = 1;
  78.     endif
  79.  
  80.     Ajj = kron (s (j:j1, j:j1), s) - eye (blksiz*n);
  81.  
  82.     rhs = reshape (b (:, j:j1), blksiz*n, 1);
  83.  
  84.     if (j1 < n)
  85.       rhs2 = s*(x (:, (j1+1):n) * s (j:j1, (j1+1):n)');
  86.       rhs = rhs + reshape (rhs2, blksiz*n, 1);
  87.     endif
  88.  
  89.     v = - Ajj\rhs;
  90.     x (:, j) = v (1:n);
  91.  
  92.     if(blksiz == 2)
  93.       x (:, j1) = v ((n+1):blksiz*n);
  94.     endif
  95.  
  96.     j = j - 1;
  97.  
  98.   endwhile
  99.  
  100.   ## Back-transform to original coordinates.
  101.  
  102.   x = u*x*u';
  103.  
  104. endfunction
  105.